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I. INTRODUCTION 

In the Standard Model (SM), CP violation in the 
quark sector is entirely described by a single weak phase 
in the CKM quark-mixing matrix. Studies of time- 
dependent CP violation in P>° —> (cc)K° decay [l[ have 
yielded precise measurements @, Q of sin 2(3, where (3 = 
arg[— {V* b V cc i) / {V t * b V t d)] and Vy are the elements of the 
CKM matrix. Measurements of time-dependent CP vio- 
lation in b — > qqs (q = u,d,s) decays offer an alternative 
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method for measuring (3. Such decays are dominated 
by b — > s loop diagrams, and therefore are sensitive to 
possible new physics (NP) contributions appearing in the 
loops of these diagrams. As a result, the effective (3 (/3 e ff) 
measured in such decays could differ from the (3 measured 
in B° — > (cc)K°. Deviations of (3 e g from (3 are also pos- 
sible in the SM, due to additional amplitudes from b — > u 
tree diagrams, loop diagrams containing different CKM 
factors ("it-penguins"), and final-state interactions. 

The decay mode P>° — > 4>K® is particularly suited for 
a NP search, as /3 e ff for this mode is expected to be very 
near in value to (3 in the SM, with sin 2(3 c g — sin 2/3 in the 
range (—0.01,0.04) However, the measurement of 

/3 e g is complicated due to other B° — > K + K~ decays 
that interfere with B° — > <j)K Q s . In general, K + K~K° S 
is not a CP eigenstate: the K + K~K® system is CP 
even (odd) if the K + K~ system has even (odd) angu- 
lar momentum. Thus, one must account for the (mostly 
S-wave) K + K~K® states that interfere with <f)K®. This 
can be done by measuring f3 c g using a Dalitz-plot (DP) 
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analysis of B° -> A+A~A°. A further benefit of a DP 
analysis is that it allows both sin2/3 e fj and cos2/3 g to 
be determined, through the interference of odd and even 
partial waves, which eliminates a trigonometric ambigu- 
ity between /3 e s and 90° — j3 e g. 

The related decay mode B + — > 4>R + is another inter- 
esting channel in which to search for NP. This decay is 
also dominated by a b — > s penguin amplitude, and its 
direct CP asymmetry, Acp, is predicted to be small in 
the SM, (0.0-4.7)% 0,0, so a significant deviation from 
zero could be a signal of NP. 

In addition to measuring /3 c g in B° —> <pK®, it is possi- 
ble to measure it for the other resonant and nonresonant 
B° — > K + K~K® decays. However, these decays may 
contain a mixture of even and odd partial-waves, so the 
final state is not guaranteed to be a CP eigenstate, thus 
posing a challenge to a measurement of j3 c g. A DP analy- 
sis can reveal which partial waves are present, thus elim- 
inating a source of systematic uncertainty affecting the 
extraction of j3 c g, without having to rely on theoretical 
predictions. 

Previous analyses of B+ -> R+R~R+ @, Q and 
B" -> K + K~K° [5l dH have revealed a complex DP 
structure that is poorly understood. Both modes exhibit 
a large peak around m{K + K~) ~ 1500MeV/c 2 , which 
has been dubbed the /x(1500). Both BaBAR and Belle 
have modeled it as a scalar resonance. The recent DP 
analysis of B° -> A° A° A° by BABAR Q2] does not yield 
evidence for this resonance. It is important to clarify the 
properties of the fx (1500) with a larger data sample, and 
in particular to determine its spin, as that affects the /3 e g 
measurement in B° — > R + R~ K®. 

An additional feature seen in B° — > K + K~K® and 
B + — > R + R~R + decays is a large broad "nonresonant" 
(NR) contribution. Previous analyses have found that 
a uniform-phase-space model is insufficient to describe 
the NR term, and have instead parameterized it with 
an empirical model. The NR term has been taken to 
be purely A+A~ S-wave in B+ -> R+R-R+ d ©], 
while smaller R + R® and A~A° S-wave terms have been 
seen in B° — > R + R~R® [To|, El, which correspond ef- 
fectively to higher-order R + R~ partial waves. Because 
the NR contribution dominates much of the available 
phase space, it is crucial to study its angular distribution 
if one wishes to accurately measure /3 e s over the entire 
B° -> R+R-R° DP. 

Because of the importance of understanding the DP 
structure in B° —> R + R~R° S) we study the related 
modes B+ -> R+R-R+ and B+ -> R Q S R° S R+ along 
with B° -> R+R-R° s . The mode B+ -> R+R-R+ is 
valuable because it has the most signal events by far of 
any B — > RRR mode. Far fewer events are expected 
in B + — > R°R°R + , but its DP has a simplified spin- 
structure due to the fact that the two K° mesons in the 
final state are forbidden (by Bose-Einstein statistics) to 
be in an odd angular momentum configuration. This im- 
plies that the /x(1500) can decay to A"A° only if it has 
even spin, and it also ensures that the nonresonant com- 



ponent in B + — > R®R®R + does not contain any A°A" 
P-wave contribution. 

In this paper we report the results of DP analyses of 
B+ -> R+R-R+ and B+ -> A°A°A+ , and a time- 
dependent DP analysis of B° -> K+K~K°. In Sec. HI 
we introduce the formalism used for the DP amplitude 
analyses. In Sec. IIII1 we briefly describe the BABAR de- 
tector and datasets used, and Sec. [IV] describes the event 
selection and backgrounds. Section[V]describes the maxi- 
mum likelihood (ML) fit parameterization and implemen- 
tation. In Sec. IVI1 we present studies of the DP struc- 
ture in the three modes, which enable us to determine 
the nominal DP models. In Sec. I VIII we then present 
the final fit results including measurements of CP vio- 
lation. We discuss systematic uncertainties in Sec. IVIIII 
and summarize our results in Sec. IIXI 



II. DECAY MODEL FORMALISM 

Taking advantage of the interference pattern in the 
DP, we measure the magnitudes and phases of the differ- 
ent resonant decay modes using an unbinned maximum- 
likelihood fit. 

We consider the decay of a B meson with four- 
momentum ps into the three daughters K\, R2, and A'3, 
with corresponding four- momenta pi, P2, and p^. The 
squares of the invariant masses are given by Sij = mfj — 

We will use the following convention for the R indices: 

• For B ± -> R ± R^R± l R 1 = R ± , R 2 = A T , and 
A3 = R^. The indices for the two like-sign kaons 
are defined such that s 12 < ^23- 

• For B ± -> A°A°A ± , R 1 = K°, A 2 = K°, and 
A3 = R^. The indices for the two A° are defined 
such that S13 < S23. 

• For B^ -> R+R-R", K x = A+, A 2 = K~, and 
A 3 = R% 

The Sij obey the relation 

S12 + S13 + s-23 = m% + m 2 Ki + m 2 K2 + m 2 K3 ■ (1) 

The DP distribution of the B ± decays is given by 

dT 1 1 (-) 2 

ds 12 ds23 = (2nf 32m 3 B+ 1 1 ' ( } 

where A (A) is the Lorentz-invariant amplitude of the 
B + (B~) three-body decay, and is a function of s\ 2 and 

S23- 

For B° — > A + A^ A°, the time-dependence of the de- 
cay rate is a function of DP location. With At = t sig —t tag 
defined as the proper time interval between the decay 
of the fully reconstructed B° -> R+R-R° s (B^ ) and 
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that of the other meson (B° ) from the T(4S), the time- 
dependent decay rate over the DP is given by 



With this definition, 



\At\/r B o 



dsi2ds23dAt 



(2tt) 3 32ml, 



4r F 



\A\ 2 + \Af 



2J^ 



E 

J<As 



(9) 



- Q (1 - 2w) (|.4| 2 - |^| 2 ) cos Am d At 
+Q (1 - 2w) 21m [e~ 2lP AA*] sin Am d Ai 



where r B o is the neutral B meson lifetime and Amj is 
the B°-B° mixing frequency. A (A) is the amplitude of 
the B s ° ig (B° sig ) decay and Q = +1(-1) when the 5° ag is 

identified as a £?° (-B ). The parameter w is the fraction 
of events in which the -B° ag is tagged with the incorrect 
flavor. 

We describe the distribution of signal events in the 
DP using an isobar approximation, which models the to- 
tal amplitude as a coherent sum of amplitudes from N 
individual decay channels ( "isobars" ) : 



'-r 1 



N 



(4) 



where 



Aj = a J F J (s 1 2,s 23 ) , 
A' = ajFj(si2,S23) 



(5) 



The i*} are DP-dependent dynamical amplitudes de- 
scribed below, and aj are complex coefficients describing 
the relative magnitude and phase of the different decay 
channels. All the weak phase dependence is contained in 
aj, and Fj contains strong dynamics only. 

The amplitudes must be symmetric under exchange of 
identical bosons, so for B + —> K + K~K + , Fj(s%2, S23) is 
replaced by Fj(sx2, S2a)+Fj(s23, S12). Similarly, in B + -t 
K°K°K + , Fj(si2,S23) is replaced by i*)(si2, S23) + 

Fj(si2,Sl 3 ). 

We parameterize the complex coefficients as 



fl 



(6) 



where Cj, bj, cf>j, and <5j are real numbers. We define the 
fit fraction (FFj) for an intermediate state as 



FFj = 



\A 3 \ 2 ) dsi 2 ds23 



JJ(\A\ 2 + \A\ 2 ) ds 12 d s 



: 23 



(7) 



Note that the sum of the fit fractions is not necessarily 
unity, due to interference between states. This interfer- 
ence can be quantified by the interference fit fractions 
FFjk, defined as 



FFik = 2 Re 



I ds^ds 



2.! 



J J(\A\ 2 + \A\ 2 ) d Sl2 dS23 



(8) 



In the B + modes, the direct CP asymmetry Acp(j) 
for a particular intermediate state is given by 



(3) AcpU) 



_ J J {\A 3 \ 2 - \Aj\ 2 ) ds 12 ds23 _ ^26, 
" J/(IA-| 2 + IA-| 2 )d S i 2 d S23 " l + b; 



(10) 



while there can also be a CP asymmetry in the interfer- 
ence between two intermediate states, which depends on 
both the b's and c)'s of the interfering states. We define 
the CP-violating phase difference as 



A(j)j = arg^j-a*) = 23 j 



(11) 



K + K Kg, we can define the direct CP 



For B° -> K + K~Kg 
asymmetry as in Eq. (|I0j) . while we can also compute the 
effective /3 for an intermediate state as 



f 



/3 offiJ - = ^arg(e 2l/3 a 3 a*) = /3 + 6 



(12) 



which quantifies the CP violation due to the interference 
between mixing and decay. 

The resonance dynamics are contained within the Fj 
terms, which are the product of the invariant mass and 
angular distributions, 

if(si2,s 23 ) - Rj(m)X L (\p*\r')X L (\q\r)Tj(L,p,q) , 

(13) 

where 

• L is the spin of the resonance. 

• m is the invariant mass of the decay products of 
the resonance. 

• Rj(m) is the resonance mass term or "lineshape" 
(e.g. Breit-Wigner). 

• p* is the momentum of the "bachelor" particle, i.e., 
the particle not belonging to the resonance, evalu- 
ated in the rest frame of the B. 

• p and q are the momenta of the bachelor particle 
and one of the resonance daughters, respectively, 
both evaluated in the rest frame of the resonance. 
For K + K~ resonances, q is assigned to the momen- 
tum of the K + , except for B~ -> K~K + K~ de- 
cays, in which case q is assigned to the momentum 
of the K~ . For KgK® resonances, it is irrelevant 
to which K° s we assign q, so we arbitrarily assign q 
to whichever K® forms the smaller angle with the 
K+. 
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• Xl are Blatt-Weisskopf angular momentum barrier 
factors [131: 



where 



L = : X L (z) = l, 
L = 1 : X z (ar) = 



1 + z 2 ' 



I 



'9 + 3zg + zfl 
9 + Zz 2 + z i 



(14) 
(15) 

(16) 



where z equals \q\r or \p*\ r', and zq is the value of 
z when the invariant mass of the pair of daughter 
particles equals the mass of the parent resonance, 
r and r' are effective meson radii. We take r' as 
zero, while r is taken to be 4 ± 2.5 ( GeV/c) -1 for 
each resonance. 



• Tj(L,p,q) are the Zemach tensors \U 
scribe the angular distributions: 



L = : Tj = 1 , 
L = l : T 3 =Ap-q, 

L = 2 : T j =^[3(p-q] 



(\p\W\? 



which de- 

(17) 
(18) 

] . (19) 



The helicity angle of a resonance is defined as the an- 
gle between p and q, measured in the rest frame of the 
resonance. For a K%K 2 resonance, the helicity angle will 
be called 9 3 , and is the angle between K 3 and K\. In 
B° — > K + K~K®, because q is defined as the K + mo- 
mentum for both B° and B° decays, there is a sign flip 
between B° and B° amplitudes for odd-L K + K~ reso- 
nances: 

F j (s 12 ,s 23 )=F j (s 1 2,s 13 ) = (-l) L ^( S i 2 , S23 ). (20) 

In contrast, for B+ -> K+R-R+ and B+ -> K°K°K+, 
Fj(si2,S23) always equals Fj(si2, 823). 

For most resonances in this analysis the Rj are taken 
to be relativistic Breit-Wigner (RBW) [15[ lineshapes: 



R 3 (m) 



1 



(m.Q — m 2 ) — imoT(m) 



(21) 



where m is the nominal mass of the resonance and T(m) 
is the mass-dependent width. In the general case of a 
spin-L resonance, the latter can be expressed as 



r(m) = r 



I h \ 2L+1 



96 



v m / 



(22) 



The symbol Tq denotes the nominal width of the reso- 
nance. The values of tuq and Tq are listed in Table [TJ 
The symbol \qo\ denotes the value of |g| when m — mo- 
For the /q(980) lineshape the Flatte form [l6[ is used. 
In this case 



Rj(m) 



(m 2 , - m 2 ) - i{g n p^(m) + gKPKK{m)) 



(23) 



p W7r (m) = yjl - 4m 2 ± /m 2 
PKK{m) = J 1 - Am 2 K /m 2 



(24) 
(25) 



Here, is the average of the and -fdJ masses, and 
g v and are coupling constants for which the values 
are given in Table [TJ 

In this paper, we test several different models to ac- 
count for NR B — > KKK decays. BaBAR's previous 
analysis Q of B+ -> K + K~K + modeled the NR decays 
with an exponential model given by 



Fnr{si 2 , S23) 



(26) 



where the symmetrization is explicit, a is a parameter 
to be determined empirically. This model consists purely 
of K + K~ S-wave decays. 

The most recently published B° — > K + K~K® analyses 
by Belle P3 and BaBAR [TTJ] both used what we will call 
an extended exponential model. This model adds K + K® 
and K~ K® exponential terms: 

A NR (s 12 , s 23 ) = a 12 e as ^ + a 13 e as ^ + a 23 e aS23 , 
A NR {s 12 , s 23 ) = a 12 e asi2 + a 13 e QS23 + a 23 e as ™ . (27) 

We also test a polynomial model, consisting of explicit 
S-wave and P-wave terms, each of which has a quadratic 
dependence on mi 2 : 

A N r(si 2 ,s 23 ) = (a so + a S ix + a S2 x 2 ) + 

(a P0 + apix + a P2 x 2 ) P\ (cos 9 3) , (28) 

where x = m\ 2 — i?, and fi is an offset that we define as 



^ = 2 \ B + \ ( mKl + mK2 + mK3 



(29) 



and P\ is the first Legendre polynomial. In this paper, 
we normalize the Pi such that 



Pe(x)P k (x)dx = 5 ik . 



(30) 



Note that in the B + — > K + K K + channel, we sym- 
metrize all terms in Eq. (|2"5)) : 

(31) 



.total = Anr(si 2 ,s 23 ) + Anr{s 2 3,si 2 ) 



This results in S-wave and P-wave terms for both the 
(K-lK 2 ) and {K 2 K 3 ) pairs. In the B+ ->■ K°K°K+ 
channel, the P-wave term is forbidden by Bose-Einstein 
symmetry. 

In Sec. I VII we present studies that allow us to deter- 
mine the nominal DP model. The components of the 
nominal model are summarized in Table [TJ Other com- 
ponents, taken into account only to estimate the system- 
atic uncertainties due to the DP model, are discussed in 
Sec. Em] 
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TABLE I: Parameters of the DP model used in the fit. Val- 
ues are given in MeV(/c 2 ) unless specified otherwise. All 
parameters are taken from Ref. [l5|], except for the /o(980) 
parameters, which are taken from Ref. [171 ]. 



Resonance 


Parameters 


Lineshape 


0(1020) 


m = 1019.455 ± 0.020 
To = 4.26 ± 0.04 


RBW 


/o(980) 


m = 965 ± 10 
= (0.165 ±0.018) GeV 2 /c 4 
g K /g n = 4.21 ± 0.33 


Flatte 


/o(1500) 


m = 1505 ± 6 
To = 109 ± 7 


RBW 


/o(1710) 


m = 1720 ± 6 
To = 135 ± 8 


RBW 


$(1525) 


m = 1525 ± 5 
To = 73±g 


RBW 


NR decays 




see text 




m = 3414.75 ± 0.31 
r = 10.3 ±0.6 


RBW 



III. THE BABAR DETECTOR AND DATA SET 



The data used in this analysis were collected with 
the BABAR detector at the PEP-II asymmetric energy 
e+e~ storage rings. The B° -t K + K~K° and B + -t 
K®K®K + modes use an integrated luminosity of 429 
ftT 1 or (471 ± 3) X 10 6 BB pairs collected at the Y{AS) 
resonance ( "on- resonance" ) . The B + — » K + K~K + 
mode uses 426 ftT 1 or (467 ± 5) x 10 6 BB pairs col- 
lected on-resonance. We also use approximately 44fb _1 
collected 40 MeV below the Y(4S) ("off-resonance") to 
study backgrounds. 

A detailed description of the BABAR detector is given 
in Ref. [l8j]. Charged-particle trajectories are mea- 
sured with a five-layer, double-sided silicon vertex tracker 
(SVT) and a 40-layer drift chamber (DCH), both op- 
erating inside a 1.5-T magnetic field. Charged-particle 
identification (PID) is achieved by combining informa- 
tion from a ring-imaging Cherenkov device and ioniza- 
tion energy loss (dE/dx) measurements from the DCH 
and SVT. Photons are detected and their energies mea- 
sured in a CsI(Tl) electromagnetic calorimeter inside the 
magnet coil. Muon candidates are identified in the in- 
strumented flux return of the solenoid. 

We use GEANT4-based [19| software to simulate the 
detector response and account for the varying beam and 
environmental conditions. Using this software, we gen- 
erate signal and background Monte Carlo (MC) event 
samples in order to estimate the efficiencies and expected 
backgrounds. 



IV. EVENT SELECTION AND BACKGROUNDS 

A. B+ -¥ K + K~K+ 

The B + — > K + K~K + candidates are reconstructed 
from three charged tracks that are each consistent with 
a kaon hypothesis. The PID requirement is about 85% 
efficient for kaons, with a pion misidentification rate of 
around 2%. The tracks are required to form a good- 
quality vertex. Also, the total energy in the event must 
be less than 20 GeV. 

Most backgrounds arise from random track combina- 
tions in e + e~ — > qq (q = u, d, s, c) events (hereafter 
referred to as continuum events). These backgrounds 
peak at cos 9t = ±1, where 0t is the angle in the 
e + e~ center-of-mass (CM) frame between the thrust 
axis of the B-candidate decay products and the thrust 
axis of the rest of the event. To reduce these back- 
grounds, we require |cos#t| < 0.95. Additional con- 
tinuum suppression is achieved by using a neural net- 
work (NN) classifier with five input variables: |cos#t|, 
| cos 6b I) |Ai/<7At|j £2/^0, an d the output of a B-flavor 
tagging algorithm. Here, 9 b is the angle in the e + e~ 
CM frame between the £?-candidate momentum and the 
beam axis, At is the difference between the decay times of 
the B + and B~ candidates with <7At its uncertainty, and 
Ck = J2j \Pj\Pk(cos9j). The sum includes every track 
and neutral cluster not used to form the B candidate, 
and 9j is the angle in the e + e~ CM frame between the 
momentum pj and the B-candidate thrust axis. Pk is the 
k th Legendre polynomial. The NN is trained on signal 
MC events and off-resonance data. We place a require- 
ment on the NN output that removes 65% of continuum 
events while removing only 6% of signal events. 

Further discrimination is achieved with the energy- 
substituted mass toes = \J (s/2 ± pi ■ ps) 2 /Ef — p B 
and energy difference AE = E B — \y/s, where (Eb,Pb) 
and (Ei,pi) are the four- vectors of the B candidate and 
the initial electron-positron system measured in the labo- 
ratory frame, respectively. The asterisk denotes the e + e~ 
CM frame, and s is the invariant mass squared of the 
electron-positron system. Signal events peak at the B 
mass (« 5.279 GeV/c 2 ) for ttt-es, and at zero for AE. We 
require 5.27 < m ES < 5.29GeV/c 2 and \AE\ < 0.1 GeV. 
An Toes sideband region with toes < 5.27GeV/c 2 is used 
for background characterization. After the calculation 
of to E s and AE, we refit each B candidate with the in- 
variant mass of the candidate constrained to agree with 
the nominal B mass [l5l |. in order to improve the reso- 
lution on the DP position and to ensure that Eq. ([1} is 
satisfied. About 8% of signal events have multiple B can- 
didates that pass the selection criteria. If an event has 
multiple B candidates, we select the one with the best 
vertex x 2 . To avoid having events that have candidates 
in both the Toes sideband and in the signal region, the 
best-candidate selection is performed prior to the toes 
and AE selection. The overall selection efficiency for 
B+ ->• K+K-K+ is 33%. 
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We use MC simulation to study backgrounds from 
B decays {BB background). In this paper, we treat 
B —> KKK decays containing intermediate charm de- 
cays as background, except for B — > XcoK (xco — > KK), 
which we treat as signal. Most of the BB backgrounds 
come from B — > D^X decays. We study 20 of the 
most prominent B + B~ background modes using simu- 
lated exclusive samples, and split these modes into six 
classes, summarized in Table |TTJ These classes have dis- 
tinct kinematic distributions, and so will be handled sep- 
arately in the ML fit, as described in Sec.fV] Class 1 con- 
tains various charmless B + decays, the largest of which 
is B + — > K + K~tt + . Class 2 includes a number of decays 
containing D° — > K + K~ in the decay chain. Class 3 in- 
cludes various decays containing D° — > K + tt~. Class 4 
consists of B+ ->■ D°K+ (D° K+ir') decays. We also 
include classes for signal-like B + — > K + K~K + decays 
coming from B + — !> D°K + (class 5) and B + J/ifjK + 
(class 6). These decays have the same Toes and AE dis- 
tributions as signal, but can be distinguished from charm- 
less signal by their location on the DP. We include a sev- 
enth BB background class, which contains the remaining 
inclusive B + B~ and B a B° decays. 



B. B + ->■ K S K° S K + 

The B + — > K®K®K + candidates are reconstructed by 
combining a charged track with two Kg — > 7r + 7r~ candi- 
dates. The charged track is required to satisfy a kaon- 
PID requirement that is about 95% efficient for kaons, 
with a pion misidentification rate of around 4%. The if' 
candidates are each required to have a mass within 12 
MeV/c 2 of the nominal K s mass and a lifetime signifi- 
cance exceeding 3 standard deviations. We also require 
that cosaK s > 0.999, where uk s is the angle between 
the momentum vector of the K ° candidate and the vector 
connecting the decay vertices of the B + and Kg candi- 
dates in the laboratory frame. The total energy in the 
event must be less than 20 GeV. 

To reduce continuum backgrounds, we require 
| cos 9t I < 0.9. We also use the same NN as for 
B + — > K + K~ K + , and place a requirement on the NN 
output that removes 49% of continuum events while 
removing 4% of signal events. Finally, the B candi- 
dates are required to satisfy 5.26 < wes < 5.29GeV/c 2 
and |A.E| < 0.1 GeV. An toes sideband region with 
toes < 5.26 GeV/c 2 is used for background characteri- 
zation. After the calculation of toes and AE, the B can- 
didates are refitted with a B mass constraint. The over- 
all selection efficiency for B + — > KgKgK + (with both 
K° s -> 7T+7T-) is 27%. 

About 2% of signal events have multiple B candidates 
that pass the selection criteria. In such cases, we choose 
the B candidate whose Kg candidates have invariant 
masses closest to the nominal Kg mass. Because there 
can be multiple B candidates that share one or more of 
the same kaon candidates, multiple B candidates may 



still remain after this step. In this case, we select the 
B candidate whose K + candidate has PID information 
most consistent with the kaon hypothesis. If multiple 
B candidates still remain, we select the one with the 
best vertex \ 2 ■ The best candidate selection is performed 
prior to the Toes and AE selection. 

BB backgrounds are studied with MC events. We 
study 10 of the most prominent background decay modes 
using simulated exclusive samples, and group them into 
three classes, summarized in Table IIII1 Class 1 con- 
tains B+ D°n + {D° -> K° S K°) and B+ -> 
K° K* + (K* + K°tv + ) decays. Class 2 contains vari- 
ous B + B~ and B°B° decays, dominated by the charmless 
decays B° -> K° S K° S K° S and B° K^ + K-R°. Signal- 
like B+ -> K° S K° S K + decays coming from B+ -> D°K + 
make up class 3. The remaining BB backgrounds are 
grouped into a fourth class. 

C. B a ^K + K-K° s 

B° — > K + K~K® candidates are reconstructed by com- 
bining two charged tracks with a Kg candidate. The 
charged tracks are required to be consistent with a kaon 
hypothesis. For most events, we apply tight kaon-PID re- 
quirements that are about 90% efficient for kaons with a 
pion misidentification rate of around 1.5%. Looser PID 
requirements (~ 95% efficient, ~ 6% pion misidentifi- 
cation) are applied in the TO12 < 1.1 GeV/ c 2 region, to 
increase the signal efficiency for B° — » (f>Kg. Kg can- 
didates are reconstructed in both the Kg — > tt + it~ and 
Kg —> ir°ir° decay modes. Kg — > it + it~ candidates are 
required to have a mass within 20 MeV/c of the nominal 
K® mass, while Kg — > 7r°7r° candidates are required to 
have a mass m^o^a in the range (m K o — 20MeV/c 2 ) < 
m^o^o < (to^o + 30 MeV/c 2 ) , where m K o is the nominal 
Kg mass. Both Kg — > tt + it~ and Kg — > 7r°7r° candidates 
are required to have a lifetime significance of at least 3 
standard deviations, and to satisfy cos a.K s > 0.999. The 
7r candidates are formed from two photon candidates, 
with each photon required to have a laboratory energy 
greater than 50 MeV and a transverse shower profile con- 
sistent with an electromagnetic shower. 

We reduce continuum backgrounds by requiring 
| cos 9t I < 0.9. In addition, we use a NN containing the 
variables | cos6*t|, | cos6>b|, and £2/ Co- Since we are per- 
forming a time-dependent analysis of B° — > K + K~ Kg, 
we omit \At/<7At \ from the NN in order not to bias the fit. 
We train the NN on signal MC events and off-resonance 
data. We make a requirement on the NN output that 
removes 26% of continuum events in the Kg — > ir + 7r~ 
channel, and 24% of continuum events in the Kg — > 7r°7r° 
channel, with only a 2% loss of signal events. B can- 
didates must satisfy 5.26 < toes < 5.29GeV/c 2 and 
-0.06(-0.12) < AE < 0.06 GeV for K° s -> tt+tt" (K° 
7r°7r°). An toes sideband region with toes < 5.26 GeV/c 2 
is used for background characterization. After the cal- 
culation of toes and AE, the B candidates are refitted 
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TABLE II: Summary of the BB backgrounds in B + — > K + K K + . The "Expected yields" column gives the expected number 
of events for 467 x 10 6 BB pairs, based on MC simulation. The "Fitted yields" column gives the fitted number of events from 
the best solution of the fit on the data (see Sec. I VILA]) . 



Class Decay Expected yields Fitted yields 



1 


B- 


— > charmless 






42 ± 5 


fixed 


2 


zr 


D M0 X,D°- 


K 


f A- 


195 ±7 


170 ± 21 


3 




h -»• D { * )0 X,D°- 


-¥ A 




117 ±5 


133 ± 34 


4 




D°K+(D°- 


> A+ 


7T-) 


92 ±5 


23 ±9 


5 


ir 


h -> D°K+(D°- 


> A+ 


A") 


233 ± 13 


238 ± 22 


6 


B- 


h ->■ J/i/)K + (J/ip - 


A 




38 ±5 


45 ± 10 


7 


B~ 


r B~ /B°B° remaining 




386 ± 12 


261 ± 56 



TABLE III: Summary of the BB backgrounds in B + — >• K®K®K + . The "Expected yields" column gives the expected number 
of events for 471 x 10 6 BB pairs, based on MC simulation. In the maximum-likelihood fit on the data (Sec. IVIIB|) . the yield 
of each class is fixed to its MC expectation. 



Class Decay Expected yields 

1 B + D°tt + {D° -> A^A"), 6.1 ± 1.2 
B+ -> K° s A*+ (A*+ ->■ K%-k+) 

2 B + /B° -> charmless 23 ± 5 

3 B + ->■ D°K + (D° -> K° S K%) 8.1 ± 1.6 

4 B+B-/B°B° remaining 118 ± 6 



with a P» mass constraint. The overall selection efficiency 
for B° -> K+K-K a s is 31% for #° -> tt+tt- and 7% for 
K° -> n°TT . 

The time difference At is obtained from the measured 
distance along the beam direction between the positions 
of the B° ig and B^ decay vertices, using the boost 
,$7 = 0.56 of the e + e~ system. We require that B can- 
didates have | At | < 20 ps and an uncertainty on At less 
than 2.5 ps. To determine the flavor of B^ g we use the B 
flavor tagging algorithm of Ref. 0], which produces six 
mutually exclusive tagging categories. We also retain un- 
tagged events (about 23% of signal events) in a seventh 
category, since these events contribute to the measure- 
ments of the branching fractions, although not to the CP 
asymmetries. 

Multiple B candidates pass the selection criteria in 
about 4% of Kg —> tt + tt~ signal events and 11% of 
K®, — > 7r°7r° signal events. If an event has multiple can- 
didates, we choose the B candidate using criteria similar 
to those used for B + -> K®K®K + . The best candidate 
selection is performed prior to the tties, AE, and At se- 
lection. 

BB backgrounds are studied with MC events and 
grouped into five classes, summarized in Table HV] We 
include classes for signal-like B° — > K + K~K® decays 
coming from B° —> D~K + (class 1), D~ K + (class 2), 
D°K° (class 3), and J/tpK° s (class 4). The remaining BB 
backgrounds are grouped into a fifth class. 



V. THE MAXIMUM-LIKELIHOOD FIT 



We perform an unbinned extended maximum- 
likelihood fit [13] to measure the inclusive B — > KKK 
event yields and the resonant amplitudes and CP- 
violating parameters. The fit uses the variables tties, 
AE, NN, mi2, and 77123 to discriminate signal from back- 
ground. Events with both charges or tag flavors Q are 
simultaneously included in the fits in order to measure 
CP violation. For B° -> K + K~K°, the additional vari- 
able At; enables the determination of mixing-induced CP 
violation. 

The selected on-resonance data sample is assumed to 
consist of signal, continuum background, and B back- 
ground components. 



A. The Likelihood Function 

1. B+ -» A+A"A+ and B + -> K° S K° S K + 

The probability density function (PDF) T>% for an event 
i is the sum of the probability densities of all event com- 
ponents (signal, qq continuum background, BB back- 
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TABLE IV: Summary of the BB backgrounds in B° — > K + K Kg. The "Expected yields" column gives the expected number 
of events for 471 x 10 6 BB pairs, based on MC simulation. The "Fitted yields" column gives the fitted number of events from 
the best solution of the fit on the data (see Sec. IVH"Cl) . 

Class Decay Expected yields Fitted yields 
Kg -» 7T + 7T- Kg -» tt'V' Kg -» tt + tt' Kg -> 7r°7r° 

1 B° ->■ D~K + (D- -> if - Kg) 42 ±13 4±1 36 ± 7 3.6 ± 0.6 

2 S° — > DJ K + (D~ — > K _ Kg) 33 ±6 3±1 11 ±4 1.1 ± 0.4 

3 B° -> D°K° S (D° -> K+K") 10 ± 1 1.0 ±0.1 16 ± 5 1.9 ±0.5 

4 B° -> J/ty>Kg(J/fy> ->■ K+K") 10 ±1 1.0 ±0.1 4±4 0.5 ± 0.4 

5 B + B~/B°B° remaining 141 ±7 123 ± 6 29 ± 28 48 ± 18 



ground), namely 

7^2 — ^sig^sig,i ± ^qq — (1 QiAqq) Pqq,i 

+ £ ^ BBj 7> (l - Q^*iy) ^BBi,i ■ (32) 
i=i 

The parameters are defined in Table El 
The PDFs T^x.i have the general form 

Px,ii s i2,s 2 3,Q) ■ (33) 

This form neglects some small correlations between ob- 
servables. Biases due to correlations in the signal PDF 
are assessed using MC events passed through a GEANT4- 
based detector simulation (see Sec. IVIII]) . 
The extended likelihood function is given by 

N 

C oc e~* l[Vi , (34) 

i 

where N is the number of events entering into the fit, 

and N = N sig + N q7} ± J2j=r N BEj is the total fitted 
number of events. 

A total of 43 parameters are allowed to vary in the 
B+ -> K+R-R+ fit. They include eight yields (sig- 
nal, continuum, and six BB background yields) and 30 
parameters for the complex amplitudes ctj from Eq. §5§ 
(see Table IVlIl in Sec. IVII[) . The last five parameters are 
A q q, one parameter each for the continuum toes and AE 
PDFs, and the means of the signal toes and AE PDFs 
(see Sec. lVCp . The A BB ~. are fixed to their MC expecta- 
tions, except for classes 5 and 6, in which they are fixed 
to the world average [l5| and 0, respectively. 

A total of 41 parameters are allowed to vary in the 
B + — > K° S K®K + fit. They include two yields (signal and 
continuum) and 16 parameters for the complex ampli- 
tudes a,j (see Table US in Sec. Iviip . The last 23 param- 
eters are A q q, one parameter each for the shapes of the 
continuum Toes and AE PDFs, the means of the signal 
toes and AE PDFs, and 18 parameters for the contin- 
uum NN PDFs (nine hot and nine f^; see Sec. lVCp . The 



A B g. are fixed to their MC expectations, except for class 
3, which is fixed to the world average fl5j |. 

2. B° -> K + K~K% 

For this decay we use a similar unbinned maximum 
likelihood fit to that described in Sec. IV A 11 but there 
are some significant differences. The components in the 
fit may be separated by the flavor and tagging category 
of the tag-side B decay. 

The probability density function Vf for an event i in 
tagging category c [2J is the sum of the probability den- 
sities of all components, namely 

V- = N sis f c V° iSii + Ny>^ + N BS f c V c BS ..(35) 

The parameters are defined in Table I VII The sig- 
nal PDF includes components for the BB background 
classes 1-4 listed in Table HVl since they lead to the same 
K+K-K° s final state. The PDFs V Xi are the products 
of PDFs for one or more variables, 

V c Xti = P X!i (m BS )P Xti (AE)P x>i (NN,s 12 ,s 23 )x 

P x /si 2 ,s 23 , At, a At, Q) , (36) 

where i is the event index. Not all the PDFs depend on 
the tagging category; the general notations P x i and V c x i 
are used for simplicity. 

The extended likelihood function evaluated for events 
in all tagging categories is given by 

7 N c 

C^Y[e-^Y[Vf, (37) 

c— 1 i 

where N c is the number of events entering into the fit in 
category c, and N c = N sig f c + N^- + N BB f c is the total 
fitted number of events in category c. 

The maximum-likelihood fit is performed simultane- 
ously over both the — » 7r + 7r~ and K® — > tt°tt q modes. 
The signal isobar model parameters are constrained to 
be equal for both modes, but otherwise the PDFs may 
differ. 

A total of 90 parameters are allowed to vary in the 
fit. They include the 18 inclusive yields (for both ifg — > 
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TABLE V: Definition of parameters in the event PDF for B + decays shown in Eq. (l32jl . The BB background classes are given 
in Tables M and |TTT] 

Parameter Definition 

iVsig total fitted B — > KKK signal yield in the data sample 

N q q fitted continuum yield 

Qi charge of the B candidate, +1 or —1 

A q q charge asymmetry in continuum events 

^ciass number of SB-related background classes considered in the fit 

N B gj fitted or fixed yield in BB background class j 

AbEi charge asymmetry in BB background class j 



TABLE VI: Definition of parameters in the event PDF for B° decays shown in Eq. I)35[l. The BB background classes are 
described in Table ITVl 

Parameter Definition 

iV s ig total fitted B — > KKK signal yield in the data sample, including BB background classes 1-4 

f c fraction of events that are tagged in category c, with ^ c f° = 1 

This fraction is assumed to be the same for signal and BB background events 
Nq-q fitted continuum yield in tagging category c 

Qi tag flavor of the event, defined to be +1 for a Stag and —1 for a St ag 

iV B g fitted yield in BB background class 5 



7r + 7r~ and K® — > 7r°7r°, there are nine yields: signal, BB, 
and seven continuum yields, one per tagging category). 
We also allow 32 parameters for the complex amplitudes 
cij to vary (22 are shown in Table IXl} six are b and S pa- 
rameters corresponding to the parameters in Table IXIII1 
and four describe the background classes 1-4 in Table HVl 
which we model as non-interfering isobars). The remain- 
ing 40 parameters include 38 parameters that describe 
the continuum PDF shapes (one AE parameter and 18 
NN parameters, for both K° -» 7r+7r~ and K° s — > 7r°7r°), 
as well as the means of the signal tties and AE PDFs for 
K° 7T+7T- only (see Sec. IV Cj) . 

B. The Dalitz Plot and At PDFs 

For B+ ->■ K+K-K+ and B+ -> K°K°K+, the signal 
DP PDFs are given by 

-P s i g (si2, S23, Q) = dT(s 12 , s 23 , Q)e(s 12 , s 23 ) , (38) 

where dT is defined in Eq. ^ , and e is the DP-dependent 
selection efficiency, determined from MC simulation. We 
assume equal efficiencies for B + and B~ events, and con- 
sider a possible asymmetry as a systematic uncertainty. 

For B° K + K~K°, the time- and DP-dependent 
signal PDF is given by 

-P s c ig (si2,S23, At, a At , Q) = 
dr{s 12 ,s 23 ,At,Q)e{s 12 , S23 ) ® K(At,a At ) , (39) 

where dT is defined in Eq. ([3]) and the At resolution func- 
tion is a sum of three Gaussian distributions. The pa- 



rameters of the At resolution function and the tagging- 
category-dependent mistag rate are determined by a fit 
to fully reconstructed B° decays 0. 

To account for finite resolution on DP location, the 
signal PDFs are convolved with a 2 x 2-dimensional res- 
olution function 

^2,^3,42,43), (40) 

which represents the probability for an event with true 
DP coordinates (s^,^^) to be reconstructed with coor- 
dinates (si2,S23)- K obeys the unitarity condition 

J fn(s r 12 ,s^ 3 ,s\ 2 ,s t 23 )ds r 12 ds r 23 = lWs\ 2 ,s t 23 . (41) 

The 1Z function is obtained from MC simulation. 

For B+ ->■ K+K-K+ and B+ -> K"K"K + , the BB 
background DP PDFs are histograms obtained from MC 
samples. The histograms have variable bin sizes calcu- 
lated using an adaptive binning method to ensure that 
fine binning is used where the DP distributions have nar- 
row structures. 

For B° K+K-K° s , the generic BB background DP 
PDFs are likewise histograms obtained from MC samples. 
The background classes 1-4 given in Table IIV1 however, 
are modeled as non-interfering isobars, and so their DP- 
and time-dependence is included in Eq. (|39[) . 

The DP PDFs for continuum events are described by 
histograms similar to those for BB backgrounds. The 
PDFs are modeled with data taken from mES sidebands, 
with a correction applied for BB backgrounds present in 
the sidebands. 
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For B° -> K+K-R", the At distribution of the con- 
tinuum events is modeled with the sum of an exponential 
decay and prompt component, convolved with a double- 
Gaussian resolution function. The parameters are taken 
from a fit to data in the toes sideband. The At distri- 
bution of the generic BB backgrounds is modeled in the 
same way, but the parameters are taken from a fit to MC 
samples. In the nominal fit, we assume zero CP violation 
in the backgrounds, but as a systematic we include CP 
violation in the BB exponential decay component. 



lenge in these fits, because some of the resonance am- 
plitudes vary rapidly as functions of the DP. The nor- 
malization of these PDFs is performed using a numeri- 
cal 2-dimensional integration algorithm that makes use of 
adaptive binning [23] . The speed of this algorithm allows 
the masses and decay widths of resonances to be varied 
in the fit. The normalization of the DP PDFs is recalcu- 
lated at each step in the fit for which these parameters 
are varied. 



C. PDFs of Other Fit Variables 

The toes an d AE distributions of signal events are 
described by modified Gaussians of the form 



P{x) oc exp 



(as - x ) 2 



2cr| + a±(x - 



x y 



(42) 



where er + and a+ are used when x > xq, and er_ and a_ 
when x < xq. Most parameters are taken from fits to 
signal MC events. The means xq are allowed to vary in 
the nominal fits to data, except for the B° — > K + K~ K®, 
K° -> n°TT° channel. 

The toes distributions for continuum events are mod- 
eled with a threshold function 2l|, while the AE distri- 
butions are modeled with first-order polynomials. The 
toes and AE shape parameters are allowed to vary in 
the nominal fits. 

A variety of PDFs are used to describe the toes and 
AE distributions of the various BB background cate- 
gories. The PDF shapes for each category are taken from 
MC simulation. Those BB backgrounds that have the 
same true final state as signal events are modeled with 
the same Toes, AE, and NN PDFs as signal events. 

The output of the NN does not have an easily pa- 
rameterized shape, so we split the distribution into ten 
bins, with the bin size chosen so that approximately 
equal numbers of signal events are expected in each bin; 
continuum events peak at larger values of the bin num- 
ber. The binned NN is then easily described using his- 
togram PDFs. The PDFs for signal and BB background 
events are taken from fits to MC events. In the case 
of B + — > K + K~ K + , due to the large number of signal 
events, we obtain the histogram bin heights for signal 
from a separate fit to data, and then fix these parame- 
ters in the nominal fit. For continuum events, the NN 
output is correlated with the distance from the center 
of the DP. To account for this correlation, the contin- 
uum NN PDF is given by a histogram with bin heights 
hi equal to hoi + giA^p- Here, Ap>p is the smallest of 

(toi2,TO 2 3,TOi 3 ). 



D. Fitting Method 

The ML fits are performed with MINUIT |22;]. Proper 
normalization of the DP PDFs poses a technical chal- 



VI. DETERMINATION OF DALITZ MODEL 

The Dalitz plots for the three B — > KKK modes are 
shown in Fig. [TJ Before fitting Acp parameters, we first 
decide which resonances and NR terms to include in the 
DP model for each of the B — > KKK modes. Because 
the B + — > K + K~K + mode has the largest number of 
events, we primarily use it to guide our decision-making, 
but the other modes are useful as well. The studies in 
this section are performed in a "CP-blind" fashion, which 
means that we constrain the CP-violating parameters of 
the signal and background components to zero, except in 
the case of B° — > K + K~ K®, where we constrain /3 e g to 
(3 for all isobars. 

One important goal is to understand the nature of 
the so-called /x(1500) resonance seen in several previous 
analyses. Both BABAR d,[ll[ and Belle [|,[l3 have mod- 
eled this resonance as a scalar particle, but while BABAR 
has found its mass and width to be inconsistent with 
any established resonance, Belle has found a mass and 
width consistent with the /o(1500). There is also confu- 
sion surrounding the branching fraction to fx(l500)K. 
Belle's B+ ->■ K+K-K+ and B° K+K~K° analyses 
both find multiple solutions for the fit fraction for fx- 
Some solutions favor a small fit fraction, less than 10%, 
while others favor a large fit fraction, greater than 50%. 
BABAR obtained a small fit fraction in B° -> K + K~K° S , 
but a large fit fraction in B + — > K + K~ K + . A large, 
broad structure around m K + K - = 1500MeV/c 2 is also 
seen by BABAR in B+ K+K--K+ [H] but not in B + H> 
K°K° s n + H|. BABAR's analysis of B° K° S K° S K° S [H 
does not provide evidence for the /x(1500). 



B 



K + K~ 



We initially perform a fit to B + — > K + K~K + using 
the same DP model as BaBAR's previous analysis of this 
mode, which includes the resonances (/>(1020), /o(980), 
/x(1500), /o(1710), and Xco, and an exponential NR 
model [Eq. (|26p]. We allow the NR parameter a, as well 
as the mass and width of the 7^(1500), to vary in the fit. 
The /x(1500) is taken to have a spin of zero. We refer to 
this hereafter as B+ K+K-K+ Model A. We find fit 
parameters consistent with BaBAR's previous analysis. 

To see how well the fit model describes the DP distri- 
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FIG. 1: Dalitz plots for B+ -> K+K~K+ (top), B+ 
K° S K° S K + (middle), and B° -> K+ K~ K° s (bottom). Points 
correspond to candidates in data that pass the full event se- 
lection, with an additional requirement that the NN output 
be 7 or less, in order to enhance the signal. 



bution, we calculate angular moments, defined as 

(P e (cos6 3 )) = J dTP e {cose 3 )dcos9 3l (43) 

where 9 3 is the helicity angle between K 3 and K\, mea- 
sured in the rest frame of K1K2, Pi is the £-th Legendre 
polynomial, and the differential decay rate dT is given in 
Eq. (T3J). Note that the angular moments are functions 
of mi2 but we suppress this dependence in our notation. 
Angular moments plotted as a function of mi2 are an ex- 
cellent tool for visualizing the agreement between the fit 
model and data, as they provide more information than 
ordinary DP projections, in particular spin information. 

If we assume that no K\K 2 partial- waves of a higher 
order than D-wave contribute, and we temporarily ignore 
the effects of symmetrization, then we can express the 
overall decay amplitude as a sum of S-wave, P-wave, and 
D-wave terms: 



<A{m\2i cos 6> 3 



A s Po(cos6 3 ) +A P e l 
+A D e l * D P 2 (cos6 3 ), 



'Pi(cos9 3 ) 

(44) 



where Ak and fa are real- valued functions of TO12, and 
we have factored out the S-wave phase. We can then 
calculate the angular moments: 



(-Pi) 
(P2) 
(P 3 ) 



V2AsAp cos (f)p + ^^J® ApAn cos (</> P 



A 2 D + V2A S A D cos 




ApAd cos (4> p - 4>d) , 



(45) 



The symmetrization of the B + — > K + K~ K + amplitude 
spoils the validity of Eq. (f4"5"j) . Nevertheless, the angular 
moments can be calculated both for signal- weighted data 
and for the fit model, providing a useful tool for checking 
how well the isobar model describes the data. In Fig. [21 
we show angular moments for data compared to the fit 
model, in the region of the DP above the (£(1020). The 
data is signal-weighted using the s Vlot [26| technique. 
The fit model histograms are made by simulating large 
numbers of events based on the fit results. In Fig. [31 we 
show the angular moments in the 0(1020) region. 

The angular moments, in particular (P2), show that 
Model A does not describe the data well in the /x(1500) 
region. If we replace the /jc(1500) with the Jo (1500) 
and the /2(1525), there is an improvement in 21n£ of 17 
units. As we will discuss shortly, this replacement is also 
motivated by a peak in (P2) seen in B + — > K®K®K + . 

We also vary the NR model. The exponential NR 
model is not very flexible; it assumes no phase motion 
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FIG. 2: B + —¥ K + K~ K + angular moments in the region mi2 > 1.04GeV/c 2 , computed for signal-weighted data, compared 
to Model A (dashed line) and Model B (solid line). The signal weighting is performed using the s Vlot method. Events with 
m K + K - near the D° mass are vetoed. 



and only an S-wave term. We fit with a polynomial 
model [Eq. ([25])] instead, which contains S-wave and P- 
wave terms and allows for phase motion. There is an 
improvement in 21n£ of 233 units. However, the poly- 
nomial model has nine more degrees of freedom than the 
exponential model. We refer to this model [which re- 
places the fx (1500) with the / (1500) and the /^(1525), 
and which uses the polynomial NR model] hereafter as 
Model B for B + — > K + K~K + . We compare the angular 



moments for Model B to data in Figs. [5] and [3] Model 
B matches the data significantly better than Model A, 
especially for (Pi) and (J-^)- 



B. B + -> K° S K^K+ 

Next we examine B + — > K®K®K + , initially includ- 
ing the resonances /o(980), /jf (1500), /o(1710), and Xco- 
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FIG. 3: B + — > K + K~ K + angular moments in the region 
rai2 < 1.04GeV/c 2 , computed for signal-weighted data, com- 
pared to Model A (dashed line) and Model B (solid line). The 
signal weighting is performed using the s Vlot method. 



We take the /x(1500) mass and width from the B + — > 
K + K~K + Model A result. We also include a polyno- 
mial NR model, but without the P-wave term, which is 
forbidden. We call this Model A for B+ -> K Q S K° S K+. 

In Fig. UJ we show the angular moments for this 
model, compared to signal-weighted data. Assuming 
there are no higher-order K^K® partial waves than D- 
wave, Eq. ([35]) is valid for B + -> K° S K° S K + . However, be- 
cause odd partial waves are forbidden in this channel, the 
odd angular moments are automatically zero. The peak 
in (P 2 ) around 1.5 GeV/c 2 in Fig. |4j suggests the presence 
of a tensor resonance. We replace the /x(1500) with the 
/o(1500) and /£(1525), and call this Model B for B+ -> 
K®K®K + . Model B improves 21n£ by 37 units over 
Model A. The angular moments for Model B are shown 
in Fig. |U Neither model does a good job of describing 
(P 2 ) in the region 1.8 < m±2 < 2.5GeV/c 2 . As an alter- 
native, we use the model from BaBAR's B° — >• K S K° S K® 
analysis 0, which includes / (980), /o(1710), / 2 (2010), 
XcO, and an exponential NR model like in Eq. ex- 
cept without the second term. For this model, 2 In C is 52 
units worse than for Model B. We then add the /2(1525) 
to this model, but its 21n£ is still 19 units worse than 
Model B. Adding the / 2 (2010) or / 2 (2300) resonance to 
Model B significantly improves 2 In C and improves the 
modeling of the (P 2 ) distribution, but no evidence for 
these resonances is seen in B + — > K + K~ K + , which has 
a much higher signal yield. Therefore, we do not include 
either of these resonances in our model. We will, how- 
ever, include these resonances as part of our evaluation 
of systematic uncertainties. 



C. B° -> K + K-K S 

Lastly, we examine B° — > K + K~ K®. We initially 
fit with the same model used in BaBAR's previous anal- 
ysis, which includes the resonances 0(1020), /o(980), 
/x(1500), /o(1710), and \c0, and the extended exponen- 
tial NR model given in Eq. (f!?T]) . We take the mass and 
width of the /x(1500) from the B+ -> K+K~K+ Model 
A result. We hereafter refer to this model as Model A for 
B° -> K + K~K°. Belle's most recent B° -> K + K~K° S 
analysis uses this same model, although with a different 
mass and width for the /x(1500). 

Using Model A, we obtain fit results consistent with 
BaBAR's previous measurement. In Fig. [5l we show the 
angular moments for this model compared to data. The 
angular moments in B° —> K + K~K® are complicated 
due to the relative minus sign between B° and B° am- 
plitudes for odd-L resonances [Eq. ([20])]. To account 
for this, when computing the odd angular moments, we 
weight the events by — Q, where Q is the flavor of the 
B t ° ag . Then, Eq. @5J) is valid for B° -> K + K-R S , ex- 
cept that for the odd angular moments, the right-hand 
side must be multiplied by (1 — 2w)/((Aro t zT B o) 2 + 1), 
which is a dilution factor caused by mistagging and B°- 
B° mixing. 
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We replace the /x(1500) by the / (1500), /a (1525), 
and / (1710), and this improves 21n£ by 18 units. We 
then replace the NR model with a polynomial NR model 
containing S-wave and P-wave terms. This improves 
21n£ by an additional 13 units. We refer to this model 
as Model B for B° — > K + K~K®; its angular moments 
are shown in Fig. [5] The improvement of Model B over 
Model A is not evident by examining the angular mo- 
ments by eye, but Model B provides a considerably better 
likelihood. 



D. Conclusion 

For each of the three decay modes, Model B produces a 
better fit to the data than Model A, at the cost of more 
free parameters. Model B also eliminates the need for 
the hypothetical /x(1500) state. The NR parameteriza- 
tion used in Model B greatly improves the fit likelihood 
in B + — > K + K~ K + , and its large number of parame- 
ters make it very flexible. A benefit of this flexibility is 
that the fit results are then less dependent on the par- 
ticular choice of NR parameterization. Model B also has 
a similar form in all three modes (the only difference is 
the absence of P-wave states in B + — > K®K®K + ), aiding 
comparison of results between the modes. In addition to 
the studies already mentioned, we tested for the presence 
of the / (1370), / 2 (1270), / 2 (2010), and / 2 (2300) in each 
mode, and in B+ -> K+R-R+ and B° -> K + K-R S , 
we tested for the 0(1680). We did not find evidence for 
any of these resonances. We also tested for the follow- 
ing isospin-1 resonances: aR(1450) in each of the three 
modes, and 0^(980) and of (1450) in B+ -> K° S K° S K+ 
and B° — > K + K~K® only. We did not find evidence 
for any of these resonances. We henceforth use Model 
B as the nominal fit model for each mode, and only in- 
clude these additional resonances to evaluate systematic 
uncertainties. 



m„ (GeV/c 2 ) 



FIG. 4: B — > K S K S K angular moments computed for 
signal-weighted data, compared to Model A (dashed line) and 
Model B (solid line) . The signal weighting is performed using 
the sPlot method. Events with m K o K o near the D° mass are 
vetoed. 
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FIG. 5: B° — > K + K~ K% angular moments in the region mi2 > 1.04GeV/c 2 , computed for signal-weighted data, compared to 
Model A (dashed line) and Model B (solid line). The signal weighting is performed using the 3 Vlot method. The plots are 
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VII. RESULTS 

A. B+ -> K + K~K+ 

The maximum-likelihood fit of 12240 candidates re- 
sults in yields of 5269 ± 84 signal events, 6016 ± 91 con- 
tinuum events, and 912 ± 54 BB events, where the un- 
certainties are statistical only. 

In order to limit the number of fit parameters and im- 
prove fit stability, we constrain the Aqp and A<f> of the 
/o(1500), /2(1525), and / (1710) to be equal in the fit 
(i.e., the b and S parameters, defined in Eq. (j6|), are con- 
strained to be the same for these isobars). We also con- 
strain the Acp and A0 of the S-wave and P-wave NR 
terms to be equal. Since the Acp in B + —> ccK + decays 
is known to be small [lj|, we fix the Acp of the \c0 to 

in the fit. Only relative values of c, (f>, and A(j) are 
measurable, so as references we fix c = 1 and <f> = for 
the NR term aso and A(j> = for all NR terms. 

When the fit is repeated starting from input parame- 
ter values randomly chosen within wide ranges above and 
below the nominal values for the magnitudes and within 
the [0 — 360°] interval for the phases, we observe con- 
vergence toward two solutions with minimum values of 
the negative log likelihood function —2 ln£ that are sep- 
arated by 5.6 units. We will refer to them as Solution I 
(the global minimum) and Solution II (a local minimum). 
The two solutions have nearly identical values for most 
parameters, but differ greatly for some of the isobar pa- 
rameters. The isobar parameters for both solutions are 
given in Table IVIII The correlation matrices of the isobar 
parameters are given in Ref. [271 ]. 

Figure |6] shows distributions of wies, AE and the NN 
output. Figure [7] shows the mi2, TO23, and 77113 distribu- 
tions for signal- and background- weighted events, using 
the gPlot [20| technique. 

The fit result for Solution I is summarized in Ta- 
ble I VIIH The systematic uncertainties are described in 
Sec. I VIIH We report branching fractions for individual 
decay channels by multiplying the inclusive branching 
fraction by the fit fractions. This neglects interference 
between decay channels. The inclusive branching frac- 
tion is computed as 

N ■ 

B{B+ -> K+K-K+) = — , (46) 

where N g g is the total number of BB pairs and e is the 
average efficiency, estimated by weighting MC events by 
the measured DP distribution, \A\ 2 + \A\ 2 . We assume 
equal number of B + B~ and B°B° pairs from the T(4S). 
We find B{B+ -> K+R-R+) = (34.6 ± 0.6 ± 0.9) x 
10~ 6 , including the XcoK + channel. We find an inclu- 
sive charmless branching fraction (excluding XcoK + ) of 
B(B+ -> K+R-R+) = (33.4 ± 0.5 ± 0.9) x 10" 6 . 

Fit fraction matrices giving the values of FFjk for So- 
lutions I and II are shown in the Appendix. Solution 

1 has large destructive interference between the S-wave 



TABLE VII: Isobar parameters (defined in Eq. ©) for B + -> 
K + K~K + , Solutions I and II. The same b and S parameters 
are used for the /o(1500), /a (1525), and /o(1710), and we 
choose to quote their fitted values in the /2(1525) rows. The 
NR coefficients are defined in Eq. ()28[) . Phases are given in 
degrees. Only statistical uncertainties are given. 



Parameter 


>». ("il 1 1 L 1 nn T 
QU1U.L1U11 1 




0(1020)^ 


c 


O.lMll ± 0.0043 


0.043 ± 0.009 






177 =b 13 


—53 zb 13 




b 


n -1- O 099 
— U.UD4 zt U.UZZ 


— U.UO( zt U.UZz 




5 


11 ± 7 


— 10 zb 6 


,/o (980)7^ 


c 


1 C<A -\- O 9Q 


1.0 zt U.O 






118 d= 12 


—34 ± 11 




b 


n nzLO -1- Oztl 

U.U'rtU ZLZ U.U^il 


— U.OZ ZIZ U.ll 




5 


4.5 ± 3.3 


-12 ± 7 




c 


U.l /9 ± U.U31 


U.28 ± U.0( 







—45 ±11 


—41 ± 15 


/K1525)K± 


c 


0.00130 ± 0.00022 


0.00160 ± 0.00038 






34 ± 10 


43 ± 16 




b 


-0.07 ± 0.05 


-0.09 ± 0.05 




6 


-0.8 ± 2.8 


0.5 ± 2.6 




c 


0.254 ± 0.044 


0.32 ± 0.08 







44 ±9 


45 ± 16 






0.114 ±0.017 


0.170 ±0.038 




6 


9 ±12 


31 ± 15 




s 


-2 ±6 


-2 ±6 


NR 










b 


-0.030 ± 0.022 


-0.062 ± 0.024 


aso 


c 


1.0 (fixed) 


1.0 (fixed) 




<t> 


(fixed) 


(fixed) 


asi 


c 


2.09 ±0.38 


0.4 ± 1.2 




4> 


160 ± 14 


1 ± 162 


as2 


c 


0.33 ± 0.08 


0.45 ±0.35 







157 ± 12 


-65 ± 19 


apo 


c 


1.6 ±0.5 


2.3 ± 1.9 




4> 


7 ±20 


130 ± 25 


api 


c 


0.80 ± 0.07 


0.85 ±0.30 




4> 


-159 ±6 


-114 ± 12 


ap2 


c 


0.49 ±0.15 


0.77 ±0.38 







-110 ± 17 


-60 ± 18 



and P-wave NR decays. Solution II has a smaller /o(980) 
fit fraction and large destructive interference between the 
/o(980) and nonresonant decays. 

We also calculate an overall charmless Acp by integrat- 
ing the charmless \A\ 2 and \A\ 2 over the DP. We find the 
charmless A CP {B+ -» K+R-R+) = (-l.7±H ± 1-4)%. 
There is negligible difference between Solutions I and II 
for this quantity. 

We plot the signal- weighted TO12 distribution sepa- 
rately for B + and B~ events in Fig. [51 Solutions I and 
II yield A C p(<)Kl020)) = (12.8 ± 4.4)% and (7.4 ± 4.5)%, 
respectively, where the uncertainties are statistical only. 
We perform a likelihood scan in Acp((f> {1020)), shown in 
Fig. [9] At each scan point, many fits are performed with 
random initial parameters, and the fit with the largest 
likelihood is chosen. Thus, the likelihood scan properly 
accounts for any local minima. The Acp is found to dif- 
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m ES (GeV/c 2 ) AE (GeV) NNbin 



FIG. 6: Distributions of m ES (left), AE (center), and NN output (right) for B + K + K~K + . The NN output is shown in 
vertical log scale. 



TABLE VIII: Branching fractions (neglecting interference), CP asymmetries, and CP- violating phases (see Eq. (|lip ) for 
B + — > K + K~ K + . The B(B + — > RK + ) column gives the branching fractions to intermediate resonant states, corrected for 
secondary branching fractions obtained from Ref. [l5|]. Central values and uncertainties are obtained from Solution I. In addition 
to quoting the overall NR branching fraction, we quote the S-wave and P-wave NR branching fractions separately. 



Decay mode 


B{B+ -> K+K-K+) x FFj (10~ 6 ) 


B(B+ -> RK+) (10~ 6 ) 


A C p (96) 


A<fe (deg) 


<^(1020)A: + 


4.48 ± 0.22±o|| 


9.2 ± 0.4±°;5 


12.8 ±4.4 ±1.3 


23 ± 13ig 


/o(980)#+ 


9.4 ±1.6 ±2^8 




-8±8±4 


9±7±6 


/ o (i500)jr+ 


0.74 ±0.18 ±0.52 


17 ±4 ±12 






fi(lS2S)K+ 


0.69 ±0.16 ±0.13 


1.56 ±0.36 ±0.30 


14 ±10 ±4 


-2±6±3 


/o(1710)A-+ 


1.12 ±0.25 ±0.50 










1.12 ±0.15 ±0.06 


184 ± 25 ± 14 




-4± 13 ±2 


NR 


22.8 ±2.7 ±7.6 




6.0 ±4.4 ±1.9 


(fixed) 


NR (S-wave) 


52±?4 ± 27 








NR (P-wave) 


24±i2 ± 27 









fer from at the 2.8 standard deviation level (2.9ct if one 
uses only the statistical uncertainties). 

Solution II exhibits a very large Aqp for the /o(980)-ftr + 
channel, but in Solution I this Acp is consistent with 0. 
A likelihood scan in Acp(fo(980)) is shown in Fig. ITUI in 
which the two solutions are clearly visible. 



B. B+ -> K° s KlK + 

The maximum-likelihood fit of 3012 candidates results 
in yields of 636±28 signal events and 2234±50 continuum 
events, where the uncertainties are statistical only. The 
BB yields are fixed to the expected number of events 
(Table ITTTj) . for a total of 155 events. 

In order to limit the number of fit parameters, we con- 
strain the Aqp and Acf> of every charmless isobar to be 
equal in the fit. We fix Acp for XcoK + to 0, but leave 
the corresponding A<f> parameter free to vary in the fit. 
Recalling that only relative values of A<p are measurable, 
our choice is therefore to measure the difference between 
A(j> for the XcO and the reference A<p shared by all the 
other isobars. 

Many fits are performed with randomly chosen starting 
values for the isobar parameters. In addition to the global 
minimum, 14 other local minima are found with values of 



—2 In £ within 9 units (3cr) of the global minimum. These 
different solutions vary greatly in their isobar parameters, 
but have consistent signal yields and values of Aqp- 

Figure QT] shows the distributions of tties, AE, and 
the NN output, compared to the fit model. Figure [T2l 
shows the TO12, TI23, and m.13 distributions for signal- and 
background-weighted events, using the s Vlot technique. 
We plot the signal-weighted m\i distribution separately 
for B + and B~ events in Fig. [T3l 

The fit result for the global minimum solution is sum- 
marized in Tables IIXI and |X] The fit fraction matrix for 
the global mininum is given in the Appendix, and the 
correlation matrix of the isobar parameters is given in 
Ref. [27j. The other minima all have consistent values 
for the /2(1525) and XcO fit fractions, but wide varia- 
tions in the fit fractions for the other states are seen. 
In particular, the fit fraction of the /o(980) varies be- 
tween 69% and 152% and the fit fraction of the / (1500) 
varies between 3% and 73%. This means the branch- 
ing fractions of these states are very poorly constrained 
with the current data. However, the signal yields for 
the different solutions only vary between 636 and 640 
events. We find a total inclusive branching fraction of 
B(B+ -> K° S K° S K+) = (10.6 ± 0.5 ± 0.3) x 10~ 6 , or 
B[B+ -+ K° 3 K° S K + ) = (10.1 ± 0.5 ± 0.3) x 10~ 6 if the 
XcO is excluded. 
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FIG. 7: Distributions of m\2 = mjf+K-,iowi m 23 = m-K+K- ,hi g iu an d m i3 — m if+K+i f° r signal-weighted (left) and background- 
weighted (right) B + — > K + K~ K + candidates in data. The event weighting is performed using the s Vlot method. The fit 
model (histograms) is shown superimposed over data (points). The signal includes the signal-like BB backgrounds (classes 5 
and 6 in Table ITT]) . The four main peaks in the upper signal plot are, from left to right: the 0(1020), /o(1500)//2(1525), D° 
(background), and XcO- The horn-like peaks in the middle and lower signal plots are reflections from the 0(1020). The XcO 
is also visible around 3.4GeV/c 2 in the middle signal plot. The upper background plot has a 0(1020) peak (mainly due to 
continuum) and D° peak (mainly due to BB). 



The global minimum has values of Acp = (4 ± 5 ± 2)% 
and Acp = (-25 ± 65 ± 11)°. The A C p for the other 
minima are between 2% and 4%. A likelihood scan of 
Acp is shown in Fig. [14] From the likelihood scan, we 
determine A CP = (4±£ ± 2)%. 



B" 



K + K~K% 



The maximum-likelihood fit is performed simultane- 
ously to 5627 candidates in the — > tt + it~ channel 
and 2910 candidates in the — > 7r°7r° channel. In the 
—> n + ir~ channel, we find 1419 ± 43 signal events 
(including 68 ± 9 signal-like BB background events, cor- 
responding to categories 1-4 in Table IIVI) . We also find 
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FIG. 8: Signal-weighted mi2 distribution for B + — > K + K~ K + candidates in data, plotted separately for B + and B~ events, 
for the entire DP range (left), and the <j!>(1020)-region only (right). The event weighting is performed using the E Vlot method. 
Signal includes irreducible BB backgrounds (classes 5 and 6 in Table HT)l . 
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FIG. 9: Scan of 2Aln£, with (solid line) and without (dashed 
line) systematic uncertainties, as a function of Acp (0(1020)) 
in B+ ->■ K + K~K + . 



FIG. 10: Scan of 2Aln£, with (solid line) and with- 
out (dashed line) systematic uncertainties, as a function of 
A C p(fo(980)) in B+ ->■ K+ K~ K+ . 



4178 ± 71 continuum events and 29 ± 28 remaining BB 
events. 

In the K° s -> 7r vr channel, we find yields of 160 ± 17 
signal events (including 7 ± 1 signal-like BB background 
events), 2703 ± 55 continuum events, and 48 ± 18 remain- 
ing BB events. All uncertainties are statistical only. 

We vary three sets of ft e g and Aqp values in the fit: one 
for the 0(1020), another for the Jo (980), and a third that 
is shared by all the other charmless isobars in order to 
reduce the number of fit parameters. Note that this last 
set of isobars contains both even-spin and odd-spin (P- 
wave NR) terms. Because of the sign flip in Eq. (f20l) . the 
sin Am^At-dependent CP asymmetry (see Eq. ([3])) has 
opposite sign for the even-spin and odd-spin components. 
We fix the f3 e g of the Xco to the SM value, and we fix its 



A CP (= -C) to 0. 

We perform hundreds of fits, each one with randomly 
chosen starting values for the isobar parameters. In addi- 
tion to the global minimum, four other local minima are 
found with values of — 21n£ within 9 units of the global 
minimum. These different solutions all have consistent 
signal yields, but vary greatly for some isobar parame- 
ters. 

Figure [15] shows distributions of tues, AE, and the NN 
output for the — > mode, and Fig. [TBI shows the 

same distributions for the K® — > 7r°7r° mode. Figure [TTl 
shows the m\2, W23, and mi3 distributions for signal- and 
background-weighted events, for the K® — > it + ir~ chan- 
nel only. Figure [15] shows the At distribution and the 
time-dependent asymmetry for signal-weighted events, 
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FIG. 11: Distributions of m E s (left), AE (center), and NN output (right) for B + -> K S K°K + . The NN output is shown in 
vertical log scale. 



TABLE IX: Isobar parameters for B + -> K° S K° S K + , for the 
global minimum. The NR coefficients are defined in Eq. (I28|l . 
Phases are given in degrees. Only statistical uncertainties are 
given. 



Parameter 


Value 


^(980)^ c 


3.35 ±0.22 
31 ±9 


^(1500)^ c 


0.20 ± 0.05 
-83 ± 18 


^(1525)^ c 


0.00179 ± 0.00032 

-58 ± 12 


MmO)^ c 


0.24 ± 0.07 
-22 ± 11 


s 


0.113 ±0.017 
45 ±60 
-12 ± 32 


NR 

b 

aso c 

4> 

asi c 

<!> 

as2 c 
<P 


-0.018 ±0.023 
1.0 (fixed) 
(fixed) 
1.00 ±0.08 

129 ± 6 
0.51 ± 0.08 
-85 ±8 



both for the 0(1020) region (1.01 < m 12 < 1.03GeV/c 2 ) 
and the 0(lO2O)-excluded region. 

The CP-conserving isobar parameters for the global 
minimum solution are summarized in Table IXI1 and the 
branching fractions are given in Table IXIII Table IXIIII 
shows the values of the CP- violating observables, with 
the central values taken from the global minimum, and 
the errors taken from likelihood scans. (Note that the 
second minimum is separated from the global minimum 
by — 2 A ln£ = 3.9, so the likelihood scan is not impacted 
by the local minima at the one standard deviation level.) 
In addition to /3 e ff and Aqp, we compute the quasi- two- 
body CP-violating parameter S, defined as 

2 Im(e- 2 ^a,a*) 1 - & 2 



The fit fraction matrix for the best solution is given in 
the Appendix, and the correlation matrix of the isobar 
parameters is given in Ref. [27] . The correlation matrix 
for the CP- violating observables is given in Table IXIV1 

The other minima all have consistent values for the 
0(1020), /2(1525), P-wave NR, and XcO fit fractions, 
but there are large variations in the fit fractions for the 
other states. Specifically, the fit fraction of the /o(980) 
varies between 19% and 41%, the fit fraction of the 
/o(1500) varies between 2% and 51%, the fit fraction of 
the /o(1710) varies between 2% and 27%, and the S-wave 
NR fit fraction varies between 34% and 120%. The signal 
yields for the different solutions, however, exhibit negligi- 
ble variation. We calculate the inclusive branching frac- 
tion using only the yield in the R® — > tt + tt~ channel. We 
find B(B° K+R-R ) = (26.5 ± 0.9 ± 0.8) x 10~ 6 , or 
B(B° R+R-R ) = (25.4 ± 0.9 ± 0.8) x 10~ 6 if the 
XcO is excluded. 

Likelihood scans for each of the /3 e g and Aqp are shown 
in Figs. [rW2T1 /3 c ff (other) is different from zero with 4.3 a 
significance. We can also distinguish between /3 e s and the 
trigonometric reflection 90° — j3 e g, due to the sensitivity 
of the DP analysis to interference between S-wave and 
P-wave amplitudes. We find that /3 e ff(other) is favored 
over 90° — f3 c s (other) with 4.8 a significance. 



D. Interpretation 

The value we measure for Acp{4>R + ) is larger than 
the SM prediction, while /3 c ff(0_ftT") is in excellent agree- 
ment with the SM. We can use the measured Acp(4>R + ) 
and /3 c g((j>Rs) to put constraints on the amplitudes con- 
tributing to these decays. We assume isospin symmetry, 
so that the amplitudes for B + <j)R + and B° 4>R° S 
are the same. We also assume that this amplitude, A, 
can be written as the sum of two amplitudes, A\ and 
A2, where A\ is the dominant penguin amplitude. A2 is 
an arbitrary additional amplitude with a different weak 
phase, which could be a tree, it-penguin, or new physics 
amplitude. 
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FIG. 12: Distributions of mi2 = m K o K a , 77123 = m if+K°,hi g hi an d m i3 — m A'+K° ,iow) f° r signal-weighted (left) and background- 
weighted (right) B + — > K%K®K + candidates in data. The event weighting is performed using the sVlot method. The fit model 
(histograms) is shown superimposed over data (points). The two main peaks in the upper signal plot are the /o(1500)//2(1525) 
and Xco- 



TABLE X: Branching fractions (neglecting interference) for B + -»■ K° S K° S K + . The B(B + ->■ RK + ) column gives the branching 
fractions to intermediate resonant states, corrected for secondary branching fractions obtained from Ref. [15J]. Central values 
and uncertainties are for the global minimum only. See the text for discussion of the variations between the local minima. 



Decay mode 


B(B + -> K° S K° S K + ) x FFj (10~ 6 ) B(B + -> i?AT+) (10" 6 ) 


/o(980)7f+ 
/o(1500)/Y+ 
/ 2 (1525)/Y+ 
/o(1710)/Y+ 

NR (S-wave) 


14.7 ±2.8 ±1.8 

0.42 ±0.22 ±0.58 20 ± 10 ± 27 
0.61 ± 0.21+S o| 2.8 ±0.9+°-^ 
0.48^ ±0.11 

0.53 ±0.10 ±0.04 168 ±32 ±16 

19.8 ±3.7 ±2.5 
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FIG. 13: Signal-weighted mi2 distribution for B + — > 
K$KsK + candidates in data, plotted separately for B + and 
B~ events. The event weighting is performed using the a Vlat 
method. Signal includes irreducible BB backgrounds (class 4 
in Table ITII)) . 
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FIG. 14: Scan of 2Aln£, with (solid line) and without 
(dashed line) systematic uncertainties, as a function of Acp 
in B+ -> K%K%K+. 



Then 



A = Ai{l + re^ +0 ), 
A = Ai{l + re^~ C) ), 



(48) 



where r is the ratio \A2/Ai\, and r\ and £ are the relative 
strong and weak phases, respectively, between Ai and 
A\ . The CP asymmetries in this case are 



and 



A CP {cj ) K + ) 



1 



2r sin £ sin r\ 



1 + 2r cos C cos rj - 



(49) 



P cS (<pK° s ) = /3±^arctan(- 



2r sin £ cos r/ + r 2 sin(2£) 



2r cos C cos 77 + r 2 cos(2£) 
(50) 

Note that Acp(cj)K^) = Acp(4>K + ) under our assump- 
tions. However, since the experimental precision on 
Acp{4>K° s ) is very poor compared to Acp{(j)K + ), we 
only include the more precise Acp{4>K + ) measurement 
in our analysis. By combining the likelihood scans of 
Acp{4>K + ) and /3 e ff {4>K®), we can put constraints on r, 
if, and (. Figure |2"21 shows the resulting constraints in the 
r-£, r-rj, and r]-( planes. 

The non-zero value of Acp(<fiK + ) leads to r — be- 
ing disfavored by 2.8a, with a value of approximately 0.1 
favored for most values of £. There is little constraint 
on £ and 77, except that values of or ±180° are dis- 
favored (because Acp(4>K + ) is non-zero), and the first 
and third quadrants of the 77-C plane are favored (because 
Acp(4>K + ) is positive). 



VIII. SYSTEMATIC UNCERTAINTIES 

The systematic uncertainties for B + — > K + K~K + , 
B+ -> K° S K° S K+, and B° ->■ K+K-K° s parameters are 



TABLE XI: CP-conserving isobar parameters (defined in 
Eq. @) for B° -»■ K + K~K°, for the global minimum. The 
NR coefficients are defined in Eq. ()28[) . Phases are given in 
degrees. Only statistical uncertainties are given. 



Parameter 


Value 


<j>(1020)K° s 


c 


0.039 ± 0.005 




<P 


20 ± 19 


/o(980)tfg 


c 


2.2 ±0.5 




<p 


40 ± 16 


/o(1500)K° 


c 


0.22 ± 0.05 




d> 


17 ±16 


#(1525)JfS 


c 


0.00080 ± 0.00028 




<P 


53 ±23 


/o(1710)JfS 


c 


0.72 ±0.11 




<P 


110 ± 11 


XcoK° s 


c 


0.144 ±0.023 




<p 


-17 ± 29 


NR 






aso 


c 


1.0 (fixed) 






(fixed) 


Osi 


c 


1.25 ±0.25 






-149 ±9 


CIS2 


c 


0.58 ± 0.22 




4> 


56 ±15 


apo 


c 


1.22 ±0.22 




<p 


65 ±13 


api 


c 


0.28 ±0.18 




<t> 


-68 ± 28 


ap2 


c 


0.42 ± 0.16 






-131 ±25 
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FIG. 15: Distributions of m E s (left), AE (center), and NN output (right) for B° -»■ K + K~K° S , K° s -> tt+tt". The NN output 
is shown in vertical log scale. 




FIG. 16: Distributions of m E s (left), AE (center), and NN output (right) for B° -> K + K~K%, K° s -¥ tt'V. The signal in the 
ttibs and AE plots has been enhanced by requiring the NN output be 6 or less. The NN output is shown in vertical log scale. 



summarized in Tables IXV[ IXVI[ and IXVII1 respectively. 
For each decay mode, the systematic uncertainties are 
assessed only for the best solution. 

We vary the masses and widths of the resonances in 
the signal model by their errors as given in Table |TJ In 
addition, we vary the Blatt-Weisskopf radii of any non- 
scalar resonances, and change the Blatt-Weisskopf radius 
of the B meson from to 1.5 (GeV/c) -1 . We take the 
observed differences in any fit parameters as systematic 
uncertainties (listed in the "Lineshape" column in Ta- 
bles |xv]|xvni • 

We vary any BB background yields that are fixed in 
the nominal fit. If the BB class contains only a single 
decay mode, the yield is varied according to the uncer- 
tainty on the world average of its branching fraction. If 
the BB class contains multiple decay modes, then we 
vary its yield by 50%. The CP-asymmetries of the BB 
background classes are also varied, either by the uncer- 
tainty on the world average or by a conservative estimate. 
Systematic uncertainties are also assigned due to the lim- 
ited sizes of the BB MC samples, which affects the BB 
PDF shapes. We also vary signal and continuum back- 
ground PDF parameters that are fixed in the nominal 
fits. This includes the parameters of the At resolution 
function and the mistag rate. An additional systematic 
uncertainty is contributed by the limited size of the data 
sideband sample used to create the continuum DP PDFs. 



These systematic uncertainties are listed under "Fixed 
PDF Params" in Tables IXVllXVlIl 

Biases in the fit procedure are studied by performing 
hundreds of pseudo-experiments using MC events passed 
through a GEANT4 detector simulation. We do not cor- 
rect for any observed biases, but instead assign system- 
atic uncertainties, listed under "Fit Bias" in Tables IXVT - 

ixvm 

We also study the effect of additional resonances that 
are not included in our nominal isobar models (see 
Sec. ED- We test for the /o(1370), ag(1450), / 2 (1270), 
/ 2 (2010), and / 2 (2300) in each mode. We also test for 
the 0(1680) in B+ -> K+R-R+ and B° -> K+K'K®, 
and the a±(980) and a J (1450) in B+ -> K°K°K+ and 
B° — > K + K~ Kg. These resonances are modeled by 
RBW lineshapes, except for the Oq (980), which is mod- 
eled by a Flatte lineshape. We first fit to data includ- 
ing these additional resonances in the model. Then, us- 
ing this fit result, we generate a large number of data- 
sized simulated datasets. We then fit to these simulated 
datasets with and without the additional resonances in 
the signal model, and take the observed differences as 
a systematic uncertainty. This is listed as "Add Reso- 
nances" in Tables IXVlXVTfl In B+ -> K + K~K + , the 
addition of the /o(1370) causes Solution II to be the 
global mininum rather than Solution I, so we do not as- 
sign a systematic uncertainty for it. 
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FIG. 17: Distributions of mi2 = rn K + K - , m^z = m K-K° s i an d m i3 = m if+Kgi f° r signal-weighted (left) and background- 
weighted (right) B° — > K + K~ K° candidates in data, K° — » tt + tv~ only. The event weighting is performed using the B Vlot 
method. The fit model (histograms) is shown superimposed over data (points). The two main peaks visible in the upper signal 
plot are due to the 0(1020) and ,/o(1500)//2 (1525). The left-most peak in the middle and lower signal plots is due to D~ /D~ 
(background). The other horn-like peaks in those same plots are reflections from the 0(1020). The upper background plot has 
a 0(1020) peak (mainly due to continuum). 



Additional systematic uncertainties are listed as 
"Other" in Tables lXVllXVIII Systematic uncertainties are 
assessed for tracking efficiency, reconstruction, and 
PID. We also compute a systematic uncertainty due 
to the limited sizes of the MC samples used to calculate 
the signal efficiency as a function of DP position. We 
assign a 1% systematic uncertainty due to possible de- 
tector charge asymmetries not properly modeled in the 
detector simulation. For the CP-violating parameters in 



B° — > K + K~ K®, we assign a systematic uncertainty 
due to the interference between CKM-favored and CKM- 
suppressed tag-side B decays [28j . 



IX. SUMMARY 

We have performed amplitude analyses of the decays 
B+ -> K+K-K+ and B+ -> K° S K°K+, and a time- 



29 





Q 

B 

B 

>-> 

< 




FIG. 18: Top: The At distributions for B° K + K~ K° s {K° s tt+tt") signal events, in the 0(1020) region (1.01 < mi 2 < 
1.03GeV/c 2 ) (left) and 0(lO2O)-excluded region (right). B° (B°) tagged events are shown as closed circles (open squares). 
The fit model for B° (B°) tagged events is shown by a solid (dashed) line. The data points are signal- weighted using the 
sVlot method. Bottom: The asymmetry (N B o — Ng ) / (N B o + Ng ) as a function of At, in the 0(1020) region (left) and 
0(lO2O)-excluded region (right). The points represent signal-weighted data, and the line is the fit model. 



TABLE XII: Branching fractions (neglecting interference) for B° -> K + K~K° S . The B(B° ->■ RK°) column gives the branching 
fractions to intermediate resonant states, corrected for secondary branching fractions obtained from Ref. 15]. In addition to 
quoting the overall NR branching fraction, we quote the S-wave and P-wave NR branching fractions separately. Central values 
and uncertainties are for the global minimum only. See the text for discussion of the variations between the local minima. 



Decay mode 


B{B° K+R-K ) x FFj 


(10" 6 ) B(B° -> RK°) (10" 6 ) 


tA(1020)A"° 
/o(980)K° 
/o(1500)/f 
f^l525)K° 
/o (1710)7,° 


3.48 ± 0.28lg;?4 

7.0±l;t ± 2.4 
0.57t°o H ±0.12 
O.lSttlt ±0-16 
4.4 ±0.7 ±0.5 


7.1±0.6±Si 

13.3+.^ ±3.2 
0-29t°oil ± 0.36 


XcoK 


0.90 ±0.18 ±0.06 


148 ± 30 ± 13 


NR 


33 ± 5 ± 9 




NR (S-wave) 


30 ± 5 ± 8 




NR (P-wave) 


3.1 ±0.7 ±0.4 





dependent amplitude analysis of B° — > K + K K®, using 
a data sample of approximately 470 x 10 6 BB decays. 

For B + — > K + K~ K + , we find two solutions separated 
by 5.6 units of — 2 ln£. The favored solution has a direct 
CP asymmetry in B+ -> 0(1020)^+ of A C p = (12.8 ± 



4.4 ± 1.3)%. A likelihood scan shows that Acp differs 
from by 2.8a, including systematic uncertainties. This 
can be compared with the SM expectation of Acp = (0.0- 
4.7)%. For B° K + R-K S , we find five solutions, 
and determine f3 c g ((f>K®) = (21 ± 6 ± 2)° from a like- 
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P J<|>(1020)) (°) 
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-100 - 
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A CP (<K1020)) (%) 



FIG. 19: Scan of 2Aln£, with (solid line) and without (dashed line) systematic uncertainties, as a function of /3 e 8 (left) and 
A C p (right) for B° -> 0(lO2O)ifg. 




50 100 
(3 eff (f o (980)) (°) 




A CP (f o (980)) (%) 



FIG. 20: Scan of 2Aln£, with (solid line) and without (dashed line) systematic uncertainties, as a function of /3 e ff (left) and 
A C p (right) for 5° -> / (980)Al!. The flat region of the fi e s scan is caused by the Acp(j o(980)) going to —100% in this region, 
in which case /3 a a becomes an irrelevant parameter. 



lihood scan. Excluding the 0(1020) K° and f (980)K° 
contributions, we measure (3 c g = (20.3 ± 4.3 ± 1.2)° for 
the remaining B° —> K + K~ decays, and exclude the 
trigonometric reflection 90° — /3 c ff at 4.8<r, including sys- 
tematic uncertainties. For B + — » K®K®K + , there is 
insufficient data to fully constrain the many complex 
amplitudes in the DP model. However, from a likeli- 
hood scan we measure an overall direct CP asymmetry 
of A C p = (4±i ± 2)%. By combining the Acp^K*) 
and p e g((j)K^) results and assuming isospin symmetry, 
we place constraints on the possible SM and NP ampli- 
tudes contributing to these decays. 

We also study the DP structure of the three B — > 
KKK modes, by means of an angular-moment analy- 



sis. This includes the first ever DP analysis of B + — > 
K®K®K + . To describe the large nonresonant contribu- 
tions seen in the three B + — > K + K~K + modes, we in- 
troduce a polynomial model that includes explicit S-wave 
and P-wave terms and allows for phase motion. We con- 
clude that the hypothetical particle dubbed the /x(1500) 
is not a single scalar resonance, but instead can be de- 
scribed by the sum of the well-established resonances 
/o(1500), #(1525), and / (1710). 
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(3 (other) (°) A cp (other) (%) 



FIG. 21: Scan of 2Aln£, with (solid line) and without (dashed line) systematic uncertainties, as a function of /3 e H (left) and 
A C p (right) for B° -> K + K~K° S , excluding the 0(1020), /o(980), and XcO- 




weak phase £ strong phase T| 




weak phase £ 

FIG. 22: Constraints on the amplitude ratio I.A2/A1.I and the relative strong and weak phases, r\ and £, between ,4i and ^2, 
for B + ->■ 07f + and S° -)■ 0ifg decays. The shaded areas, from light to dark, show the la, 2<r, 3er, and 4<r allowed regions. 
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TABLE XIII: CP- violating parameters /3 c s, Acp, and S for B —¥ K + K Kg. Central values correspond to the global 
minimum. Statistical uncertainties for /3 e ff and Acp are determined from likelihood scans. 



Component 


/3 cff (dcg) A CP (= -C)(%) 


S 


0(1020)/^ 
/o(980)tf° 
Other 


21±6±2 -5 ±18 ±5 
18 ± 6 ± 4 -28 ± 24 ± 9 
20.3 ±4.3 ±1.2 -2±9±3 


0.66 ±0.17 ±0.07 
0.55 ±0.18 ±0.12 
0.65 ±0.12 ±0.03 



TABLE XIV: Statistical correlation matrix for the CP- violating parameters f3 c g and Acp for B° — > K + K K%. The matrix 
corresponds to the global minimum solution. 



/3 eff (0(1020)) 
/3eff(/o(980)) 
/3 eff (Other) 
A C p(<p(W20j) 
A C p(fo(980)) 
A C p (Other) 



/3 cff (^(1020)) /? off (/o(980)) /3 cff (Other) ^(0(1020)) Acp(/o(980)) Acp(Othc 



1.00 



0.38 
1.00 



0.15 


0.21 


-0.44 


-0.32 


0.63 


-0.10 


0.05 


-0.33 


1.00 


-0.13 


0.47 


0.14 




1.00 


-0.25 


-0.14 






1.00 


0.60 








1.00 



TABLE XV: Summary of systematic uncertainties for B + — > K + K K + parameters. Errors on phases, Aqp's, and branching 
fractions are given in degrees, percent, and units of 10~ 6 , respectively. 



Parameter 


Lineshape Fixed PDF Params Other Add Resonances Fit Bias Total 


A0(0(1O2O)) 


3 


1 





2 


2 


4 


A0(/ O (98O)) 
A0(/£(1525)) 


2 


1 





6 


1 


6 


1 








3 


1 


3 


A0(xco) 


1 


1 





1 


1 


2 


A C p (0(1020)) 


0.2 


0.2 


1.0 


0.3 


0.7 


1.3 


Act (/o(980)) 


3 


1 


1 


2 


1 


4 


A C p (/a (1525)) 


1 


1 


1 


3 


1 


4 


Acp (NR) 


1.1 


0.4 


1.0 


0.8 


0.7 


1.9 


5(0(1020)) 


0.20 


0.04 


0.11 


0.14 


0.08 


0.29 


B(/o(980)) 


1.2 


0.1 


0.3 


2.5 


0.4 


2.8 


B(/o(1500)) 
B(/a(1525)) 


0.06 


0.02 


0.02 


0.52 


0.02 


0.52 


0.05 


0.01 


0.02 


0.07 


0.10 


0.13 


B(/o(1710)) 


0.08 


0.04 


0.03 


0.49 


0.05 


0.50 




0.01 


0.01 


0.03 


0.02 


0.04 


0.06 


S(NR) 


1.0 


0.2 


0.5 


7.4 


0.3 


7.6 


B(NR (S-wave)) 


13 


2 


1 


23 


2 


27 


B(NR (P-wave)) 


10 


2 


1 


25 


3 


27 


B(Total) 


0.0 


0.2 


0.8 


0.1 


0.4 


0.9 


ZJ(Charmless) 


0.0 


0.2 


0.8 


0.1 


0.3 


0.9 



TABLE XVI: Summary of systematic uncertainties for B" 1 
are given in percent and units of 10~ 6 , respectively. 



KsK®K + parameters. Errors on Acp and branching fractions 



Parameter 


Lineshape Fixed PDF Params Other Add Resonances Fit Bias 


Total 


A C p 








1 





1 


2 


B(/o(980)) 


1.4 


0.3 


0.3 


1.0 


0.4 


1.8 


B(/o(1500)) 
B(/£(1525)) 


0.05 


0.03 


0.01 


0.57 


0.04 


0.58 


0.06 


0.02 


0.02 


0.07 


0.03 


0.10 


B(/o(1710)) 


0.06 


0.04 


0.01 


0.02 


0.08 


0.11 


B(Xco) 


0.01 


0.01 


0.01 


0.00 


0.03 


0.04 


23(NR (S-wave)) 


1.3 


0.6 


0.4 


2.0 


0.2 


2.5 


B(Total) 


0.0 


0.2 


0.2 


0.0 


0.0 


0.3 


£>(Charmless) 


0.0 


0.2 


0.2 


0.0 


0.0 


0.3 
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TABLE XVII: Summary of systematic uncertainties for B° — > K + K Kg parameters. Errors on angles, Acp's, and branching 
fractions are given in degrees, percent, and units of 10~ 6 , respectively. 



Parameter 


Lineshape Fixed PDF Params Other Add Resonances Fit Bias Total 


/3 off (0(lO2O)) 


2 


1 





2 





2 


/3 eff (/o(980)) 


1 


1 





4 





4 


/3 e fi (other) 


0.7 


0.4 


0.2 


0.8 


0.4 


1.2 


A C p (0(1020)) 


2 


2 


2 


2 


3 


5 


A C p (/o(980)) 


6 


3 


2 


5 


2 


9 


Acp (other) 


1 


1 


1 


2 


1 


3 


£(0(1020)) 


0.13 


0.05 


0.08 


0.05 


0.03 


0.18 


B(/o(980)) 


1.3 


0.3 


0.1 


2.0 


0.1 


2.4 


B(/o(1500)) 
B(/^(1525)) 


0.04 


0.02 


0.02 


0.10 


0.03 


0.12 


0.02 


0.01 


0.00 


0.15 


0.02 


0.16 


B(/o(1710)) 


0.3 


0.1 


0.1 


0.4 


0.1 


0.5 


B(Xco) 


0.02 


0.02 


0.02 


0.01 


0.04 


0.06 


B(NR(Total)) 


2 


1 


1 


8 


1 


9 


B(NR (S-wave)) 


2 


1 


1 


8 


1 


8 


B(NR (P-wave)) 


0.1 


0.2 


0.1 


0.3 


0.1 


0.4 


B(Total) 


0.0 


0.4 


0.7 


0.0 


0.1 


0.8 


^(Charmless) 


0.1 


0.4 


0.6 


0.0 


0.2 


0.8 
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APPENDIX 

We give tables of the interference fit fractions FFjk, 
denned in Eq. ©. 
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TABLE XVIII: Values of the interference fit fractions FF jk for B + ->■ K + K~K + , Solution I. The diagonal terms FFjj are the 
ordinary fit fractions FFj, which sum to 272%. The NR component is split into S-wave and P-wave parts for these calculations. 
Values are given in percent. 



0(1020) 
/o(980) 
/o(1500) 
/a (1525) 
/o(1710) 

XcO 

NR (S-wave) 
NR (P-wave) 



0(1020) /o(980) /o(1500) / 2 (1525) /o(1710) XcO NR (S-wave) NR (P-wave) 



12.9 



-0.1 
27.2 



0.0 


0.0 


0.1 


-0.0 


-7.4 


8.2 


-4.7 


-0.0 


-5.4 


-1.0 


-0.8 


-3.7 


2.1 


0.0 


2.3 


0.1 


3.1 


-0.8 




2.0 


0.1 


-0.0 


-0.0 


0.7 






3.2 


-0.1 
3.2 


-13.5 
3.3 
151.4 


4.9 
-1.8 
-155.0 
69.4 



TABLE XIX: Values of the interference fit fractions FF jk for B + — > K + K~ K + , Solution II. The diagonal terms FFjj are the 
ordinary fit fractions FFj, which sum to 174%. The NR component is split into S-wave and P-wave parts for these calculations. 
Values are given in percent. 





0(1020) / (980) /o(1500) / 2 (1525) /o(1710) 


XcO 


NR (S-wave) NR (P-wave) 


0(1020) 


12.3 -0.3 -0.1 -0.0 


-0.1 


-0.1 


-1.5 


5.1 


/o(980) 


12.5 1.5 0.1 


3.9 


0.6 


-40.6 


-10.2 


/o(1500) 
/ 2 (1525) 
/o(1710) 

XcO 


2.6 -0.0 
1.5 


2.3 
0.0 
2.5 


0.1 
-0.0 
-0.0 

3.6 


-3.5 
-0.3 
-11.6 
-1.5 


-0.0 
0.7 

-2.4 
0.5 


NR (S-wave) 








91.1 


-17.2 


NR (P-wave) 










48.2 



TABLE XX: Values of the interference fit fractions FF jk for B + -> K^K° S K + , for the global minimum. The diagonal terms 
FFjj are the ordinary fit fractions FFj, which sum to 345%. Values are given in percent. 



/o(980) 
/o(1500) 
/ 2 (1525) 
/o(1710) 

XcO 

NR (S-wave) 



/u(980) /o(1500) / 2 (1525) /o(1710) X co NR (S-wawf 



139.0 



-19.2 
4.0 



0.0 
-0.0 
5.7 



-12.4 


-1.0 


-217.0 


4.1 


0.2 


9.5 


-0.0 


-0.0 


-0.0 


4.5 


0.1 


-9.2 




5.0 


-0.0 






186.5 



TABLE XXI: Values of the interference fit fractions FF jk for B° -> K + K~K° S , for the global minimum. The diagonal terms 
FFjj are the ordinary fit fractions FFj, which sum to 188%. The NR component is split into S-wave and P-wave parts for these 
calculations. Values are given in percent. 



0(1020) 

/o(980) 

/o(1500) 

/ 2 (1525) 

/o(1710) 

XcO 

NR (S-wave) 
NR (P-wave) 



0(1020) /o(980) /o(1500) / 2 (1525) /o(1710) XcO NR (S-wave) NR (P-wave) 



13.1 



0.0 
26.3 



0.0 


0.0 


0.0 


0.0 


0.0 


0.2 


0.1 


-0.0 


14.4 


-0.7 


-81.2 


0.0 


2.1 


-0.0 


5.3 


-0.1 


-0.7 


0.0 




0.5 


-0.0 


0.0 


0.0 


0.0 






16.7 


-0.2 


-27.0 


0.0 








3.4 


1.6 


0.0 










114.5 


0.0 












11.7 



